Computationally efficient phase-field models with interface kinetics 
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We present a new phase-field model of solidification which allows efficient computations in the 
regime when interface kinetic effects dominate over capillary effects. The asymptotic analysis re- 
quired to relate the parameters in the phase-field with those of the original sharp interface model 
is straightforward, and the resultant phase-field model can be used for a wide range of material 
parameters. 
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Phase-field model techniques have become increasingly 
recognized as the tool of choice for solving moving free 
boundary (sharp-interface) problems, and in particular 
solidification processes 0, S I3> 13 ■ Following early work 
that related the phase-field models to the original sharp 
interface model in the limit of zero interface width ||5|, 
Karma and Rappel showed that calculations could be 
performed when the interface width is of the order of the 
capillary length Q , even scaling with the dendritic tip 
radius at low undercooling. This work provided a par- 
ticular relation between the parameters of the original 
sharp interface model and the phase-field model, and has 
been the starting point for accurate computations of so- 
lidification in the important regime when capillary effects 
dominate interface kinetics @. 

Recently however, there has been a growing interest 
in the opposite regime where interface kinetics are dom- 
inant. This interest is stimulated by experimental obser- 
vation of the puzzling morphological transition of the so- 
lidification front of Ni at high undercooling jj flollTlLIT^ . 

The purpose of this paper is to present a modification 
of the phase-field model of solidification so as to enable 
efficient computations in the regime when interface ki- 
netics are the dominant factor. The modification allows 
one to use an interface thickness many times larger than 
the capillary length. The methodology for performing 
the asymptotic analysis is different and much easier to 
perform systematically to high order than the asymp- 
totic matching employed in all previous analyses, and is 
capable of being used in other free boundary problems. 

Sharp interface model:- The symmetric model for the 
solidification of a pure melt from the liquid (L) phase to 
the solid phase (S) is defined by the equations: 



dnU\s - dnU\L = V/D 

u, + dok = -B{V) . 



(1) 
(2) 
(3) 



Here u = {T — Tm)c/ L is the dimensionless temperature 
and Ui is its value at the solidification front, with T be- 
ing the temperature in the liquid or solid, Tm being the 
melting temperature of a planar interface, c being the 



specific heat and L being the latent heat of fusion per 
unit volume. The curvature of the solidification front is 
given by k and the capillary length is do- D is thermal 
diffusivity (assumed here to be the same in both phases), 
and V is the normal velocity of the front. Equation Q 
expresses heat conservation at the interface, and equa- 
tion (PI is a modified Gibbs- Thomson condition, which 
is a statement of local equilibrium at the interface with 
the attachment kinetics included through the term B{V). 
Traditionally, a linear kinetic undercooling B{V) — (3V is 
used. It should be stressed, however, that the linearity of 
Ui with respect to velocity is a purely phenomenological 
assumption, and molecular dynamics simulations |l3lll4| 
suggest that there are substantial deviations from it at 
large undercooling. In this paper we are interested in ma- 
terials with large dimensionless parameter /3 = PD/dQ. 
This constant is a measure of the importance of interface 
kinetics for a given material. It takes very different val- 
ues for different materials, and for Ni it is estimated from 
molecular dynamics simulations to be as high as 90 [l5l |. 

The phase-field model:- The phase-field equations can 
generally be written in the form: 



dtu = DW^u+l-dthii;). 



(4) 
(5) 



Here ip represents the phase-field, f{ip) is a double well 
potential, g{ip) shifts the relative height of the two min- 
ima making one of the phases metastable for u 7^ 0. Note 
that we have used the subscript notation to denote differ- 
entiation: hence gj/,('0) means dg/dip. The sharp inter- 
face boundary is recovered as the locus of points where 
'0 = 0, and we are interested in the behavior of the 
phase-field equations as the phase-field interface width 
W and relaxation time r tend towards 0. In order to 
solve the desired sharp interface model, we need to as- 
certain what phase-field parameters {r, W, A, f {i^) , gii^)} 
should be used given the values of {D, do, (3 or B{V)}. 

Asymptotic analysis:- The limit of small ly is a singu- 
lar one because W multiplies a highest order derivative. 
All previous works use asymptotic matching to deal with 
the singularity. Here we will demonstrate a simpler ap- 
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proach based on the fact that the equation for u is Hnear 
so that it can be trivially solved in terms of ■0- Following 
Karma and Rappel, the analysis we give here focuses on 
the coordinate perpendicular to the interface, r; the effect 
of the transverse dimensions is incorporated by curvature 
dependent corrections, just as in their work. 

Requiring that u is finite at r — > ±00 and requiring 
that i^ir) ^ t1 sufhciently fast one obtains 



one wishes to evolve the system. For a spatially explicit 
numerical scheme At < ir {Ax/W)^ /A. The factor of A 
is included to guarantee accuracy in the presence of the 
term Xug^. Collecting terms, and using W ~ W cc Ado 
and r oc A^ we obtain 



ocAM — 
W 



(12) 



u{0 = + l-pe 



1 1 

+ -p\u^ — ■ — 

2 V p + <i 



.ip+i)vh{il;{T]))dr] (6) 



(7) 



where 



Vt 

kW, w = — < 1, 

w 



vw 

D ' 



^ = r/W, (8) 



m = sgn{p + q), and the last line defines ii. The de- 
pendence on p + q expresses the singular nature of the 
problem with respect to this parameter. From u(^) one 
can compute the outer limit Uoutif) and thus obtain 
Ui = Uout{0) = uo + 1/2 p/{p + q)- 

Despite the singularity, if we are interested only in the 
profile of u near the interface one can expand in powers 
of p and p + q and to first order obtain 

1 /■« 

U{0 = U^ + -P J (hi^Piv)) ~ 1) dTJ + O {p{p + q)) . (9) 

Using this expression and substituting it back in the 
equation for ip we get an equation for ip which can be 
solved using regular perturbation theory. In this way 
we recover the standard asymptotic result of Karma and 
Rappel: 



WiX) = A- 



r(A) = A^ 



I3D 
do 



a2X 



dl 



(10) 



where ai and 02 are constants depending on /(ip) and 
g{ip). Notice also that since ip is monotonic one can com- 
pute the distance from the interface ^ from ip and in this 
way write to order p 



1 



(11) 



where Fi is defined as the integral in equation 

From H10|) follows that with the functions /(•), g(-), 
and h{-) fixed, A is the only free parameter. 

Computational complexity:- We now examine how the 
computation time tc for the phase-field model scales with 
the free parameter A in a discretized calculation with 
adaptive mesh refinement and uniform grid elements in 
d dimensions. Clearly, tc depends on the width of the 
phase-field boundary layer W, the space resolution Ax, 
and the time step At. The inverse computation time 
scales as: {At/t){Ax/W){Ax/L)'^-^ where L'^~^ is the 
order of the surface area and t is the maximum time 



showing that the computation time is highly sensitive to 
A and also depends on the spatial resolution required by 
the shape of the interface profile of ip: smoother profiles 
are better! 

While it would be computationally efhcient to work 
with large A, doing so would introduce higher order 
terms in the curvature k and velocity V into the Gibbs- 
Thomson condition Q . In principle, it might be possible 
to fine tune g{^p), f{ip), and h{ip) to kill terms of order 
p^, pq, w^, vq, etc. However, the resultant expressions are 
very complicated, the integrals involved cannot be done 
analytically, there are many terms to consider, and even 
if successful, this would most likely introduce delicate fine 
structure into the phase-field profile which would offset 
the computational benefits. 

How can we do better? Using equation (|10|) we see 
that we require v — {pi + a2X)p ^ 1 which puts a very 
severe constraint on p if /3 is large. Computationally this 
would be more important at large undercooling when a 
thin temperature boundary layer forms around the so- 
lidification front, and correspondingly q < p so that the 
smallness of p is the limiting factor. (For example at un- 
dercooling A = 0.8 the theory predicts q/p = 1/7 at a 
steady state dendrite tip.) Therefore, a good objective 
is to modify the phase-field in a way that relaxes the 
constraint on v. 

A step in that direction was made by Bragard et al. 
[Tsf . who replaced Xu by H{Xu) in equation iJ(-) is 
computed numerically by solving the following non-linear 
eigenvalue problem with appropriate boundary condi- 
tions on ip: 



d^i} 



d^/j 



U{i')+v^~H{v)g^{i,)=Q . (13) 



dx 



With H{-) chosen in this way the non-linearities appear- 
ing in the standard phase-field model at large v are can- 
celled by the non-linearities in _ff(-). To relate the pa- 
rameters they use 



W T 



(14) 



However, this relationship is valid only in the limit of van- 
ishing p, i.e. when u is approximately constant across the 
diffuse interface - a result analogous to that of Caginalp 
for the standard phase-field. Correspondingly, in |l3| a 
value oi p close to 0.01 is used in computations. 

Since t; -I- g is no longer a small parameter it is analyt- 
ically more involved to derive corrections for finite p. To 
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compute the linear part it is enough to consider small v 
and the result is 



W 

r = \W{P + a2 — ) 



(15) 



which is the analog of Karma and Rappel's formula. 

Replacing H14|l by H15() allows much larger values of p 
to be used and is a straightforward way to make better 
use of the model proposed in Numerically we ob- 

served that for = (1 — ^^)^, the form used in [l5l| . 
the non-linearities in k and V are weak even for values 
of V of the order of 20. However, the Bragard et al. 
phase-field model, even with the improved asymptotics 
(|15|l that we derived, still does not provide the desired 
degree of computational improvement because the phase- 
field profile develops a new length scale of order W/H{v) 
which needs to be resolved numerically in order to avoid 
artifacts. In addition, H increases very rapidly with v. 

New class of models:- We propose to replace r by Tfi{il}) 
in such a way so that the effective equation for ■0 becomes 



(16) 



What are the advantages of doing this? The asymptotic 
analysis is greatly simplified because the equation for V' 
can be analyzed separately from that for u. The solution 
for is simply '(/;(^) = V'o(f) where V'o(f) is the solu- 
tion of 91-00 — fipii^o) = 0, and the relation between the 
parameters is simply given by H14() . 
Now we can rewrite (|16|) as 

rdt^p = W^V^^-f^--XuW\Vij\~X^Fi{^p)W\V^j\ (17) 

with the only problem being the presence of p in the 
evolution equation. The final trick is to express p in 
terms oi dtip: 

^ D D da^ij D IW-I ■ ^ ' 

The equation for ijj becomes 

TR^dti^ = W^'VV - UW - XuW\V^Pl (19) 



where 



1 



(20) 



For f{4>) = I (l — tJj'^Y ^'^d f^i^) = V' we have Fi{'ip) — 
V21n((-0 -I- l)/2). It follows that tr > t which means 
that the model is well behaved. The expression for tr can 
be compared with Karma and Rappel's formula which 
can be rewritten in the form t' — r + 02 A 02 is 
approximately the value of — l/2i^i(0). As it stands 
tr(—1) = 00 and points with tp = —1 cannot evolve, 
so some cutoff near ip = —1 should be introduced. Ex- 
periments show that results are insensitive to the exact 
form of the cutoff. 



The restriction of order p accuracy comes from expand- 
ing u(^) near the interface. It is possible to go to higher 
orders in p or simply use the full expression This 
would result in an implicit equation for dtip 



2D 



dtipu{'ip,p + q). 

(21) 

The function ii can be tabulated in advance and equation 
(|21|1 can be solved iteratively at each time step. 

Different approximations to eqn. (|21|l lead to different 
schemes. For example if we consider the next order term 
in ^p{p + q)F2 {ip\ {V''(C)}) ^nd up with a quadratic 
equation for dtip!. The model including p^ corrections is: 

TR = T->^{F^{^)+qF2m 
rR{dti^)o = W^W^^ - - XWu \V^\ 

a - (CtV')0 7^2D^ g^W 



dti! = (5t?A)o 



(at0)o(l + a + 2a2 



.X.22) 



In the evolution equation q = Wk ~ WV ■ n with n = 

VV'/IVV'I. 

Another application of the technique is to correct for 
terms of order pq at small undercooling when p <ti q. In 
this regime terms of order pq are not negligible compared 
to terms of order p. To achieve this it is enough to use 
(9*0)0 from above without correcting it. 

The phase-field model can be generalized to handle ar- 
bitrary interface kinetics Ui = —dok—B{V) at orderp and 
arbitrary v. The resultant phase-field model equation is 

TR{i;,u)dt^ = W^2(vV-fc(V^-n))-/v,(V') 

- B-'^ {X{u + dok))W\Wi;\ (23) 
XW^ 



tr{iP,u) = t- 



2D 



■Fi (0)6-1 iX{u + dok)){2A) 



The above recipe for improving phase-field models can 
be used also in cases when the ^ profile changes with V 
and k. For example it can be applied to the model of 
Bragard et al. by effectively replacing u with Ui yielding 

rfl,(0, u)dti^ = W^V^i; - U{iP) - i/(-AM)<7^(0), (25) 



with 



TR 



r + ^il'(-A.)A(^,i/(-A.))^. (26) 



In this case the expression for tr is more complicated 
because tp changes with v. The functions H{v) and ipiviC) 
which solve equation (|13|l and Fi can be pre-computed 
numerically leading to a very efficient numerical scheme. 

Anisotropy:- To include anisotropy we need only to re- 
place W with W{n), T with T(n), and in three dimen- 
sions, 

W^2v20 with - ^ / dFVK(n)2(V0)2 to obtain the 
Gibbs-Thomson condition 



l E [»-(■., + »■(■»] i--^V-. (27, 



i=l,2 
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FIG. 1: Interface velocity from different phase- field models 
as a function of A. /3 = 90, A = 1.2. At t = 8 x lO^di/D. 
Ax = 0.5W = O.SAdo 



where 9i and ^2 are the angles between the normal and 
the local principle directions on the interface, and R^^ 
and i?^^ — the principle curvatures. 

Numerical Experiments:- We now compare the perfor- 
mance of different phase-field models in one-dimensional 
simulations. The benchmark problem solved is u(t — 
0,x) — — A for X G (— cx3,oo) with the solid-liquid inter- 
face initially at x = 0. The interface velocity, V{t), is 
compared to that for the sharp interface model obtained 
via direct numerical integration. 

The models compared are identified as follows. ST: 
standard phase-field model Q; BR: Bragard et al. 
model with asymptotic relation H14I) : BR+: the above 
model with the improved asymptotic relation l|15|) : t^: 
the new model l(T^ : tr+p^: the new model with cor- 



rections (|22() : tb BR: the improved version of BR given 
in We used h{^) = ip and = (1 - V'^)^. 

As an example, we computed the velocity of the front 
after t = 3.5 x lO'^dl/D for /3 10, A = 1.2 and A = 15. 
The exact result is Vdo/D — 0.021. The results for mod- 
els ST, BR, TR BR were 0.012, 0.044, 0.020 respectively, 
showing that the previously existing models are inade- 
quate in this regime. 

Figure ^ compares the systematic deviations of vari- 
ous phase-field models from the sharp interface solution 
as a function of A. One clearly sees that BR model 
leads to errors linear in A. For BR+ unintended non- 
linearities in the Gibbs-Thomson condition quickly in- 
crease the error with A. In contrast tr and tr + p^ 
models yield approximately the same velocity for the en- 
tire range of A considered. Using the values for Ni cited 
in 15], D = lO^^mVsec and do = 5.56 x 10"^", the 
choice A = 1.2 corresponds to a steady state velocity 
V = {A — l)D/(/3do) = 40m/sec which is approximately 
where the experimentally observed morphological transi- 
tion occurs. 

To match the accuracy of tr model with A = 30 we 
need to take about A = 5 in BR (measuring deviations 
from the limiting phase field value). In 3D, this leads 
to about (30/5)^ ^ 200 times increase in computational 
speed as compared to the simulations in [l5l | . If we use 
the Bragard et al. model with our improved asymptotics, 
the new tr{iP) models will be about 3'^ = 27 times faster. 
The above figures are just for illustration, the precise 
computational gains will depend on the desired accuracy 
and the regime of interest. 

In conclusion, the models described here are the first 
that can systematically handle interface kinetics domi- 
nated growth in and beyond the thin-interface limit en- 
abling huge gains in computational efficiency. 
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